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SECTION  I 


INTRODUCTION 

A  hypervelocity  arc  tunnel  can  be  described^  as  a  convergent -divergent  duct 
characterized  by  (a)  a  relatively  short  plenum  chamber  into  which  high 
enthalpy,  high  pressure  gases  are  introduced  from  an  arc-heater  section 
(possibly  mixed  with  a  secondary  gas  flow)  at  subsonic  velocity,  and 
typically  at  or  near  thermochemical  equilibrium  (but  with,  perhaps,  vibra¬ 
tional  excitation),  (b)  an  area  contraction  into  the  throat  region,  and 
(c)  a  nozzle  section  through  which  the  flow  expands  to  supersonic  or  hyper¬ 
sonic  velocity,  with  concomitant  reductions  in  pressure  and  temperature, 
at  times  (depending  on  geometry)  producing  a  rapid  slowing  and/or  freeze- 
out  of  the  chemical  reactions  and  vibrational  relaxation  in  the  nozzle 
section  .  Due  to  the  relatively  short  length  of  the  plenum  section,  the 
nature  of  the  arc  heating  process,  and  the  possibly  incomplete  mixing  of 
primary  and  secondary  gas  streams,  significant  inhomogeneities  may  exist 
in  the  gas  flow  entering  the  throat  region,  both  in  terms  of  gas  composi¬ 
tion,  velocity  and  thermodynamic  state.  These  inhomogeneities  will 
convect  through  the  nozzle,  producing  uncertainties  and  possibly  spurious 
effects  in  the  interpretation  of  test  results.  The  requirement  for  the 
computer  program  which  has  been  developed  arises  from  the  need  to  better 
understand  how  these  inhomogeneities  convect  through  the  tunnel,  and  to 

quantitatively  assess  their  effects  on  the  uniformitv  of  the  flow  produced 
at  the  nozzle  exit  (or  test  section). 

Thus,  the  flow  problem  to  be  numerically  solved  consists  of  a  nonuniform 
two-dimensional  or  axisymmetric  high  enthalpy  flow  which  may  be  in  or 
near  equilibrium  in  a  plenum,  and  expands  through  an  arbitrary  duct.  For 
the  program  to  be  useful  for  enthalpies  up  to  10*4  BTU/lbm  and  pressures 
to  200  atmospheres,  a  vibrational  energy  equation  must  be  included  in  the 
system  of  equations  for  each  diatomic  or  polyatomic  species,  as  their 
vibrational  states  may  be  character ized  by  temperatures  significantly 
different  than  the  static  temperature.  In  the  axisymmetric  case,  a  swirl¬ 
ing  flow  is  often  used  to  stabilize  the  arc,  resulting  in  an  appreciable 
level  of  angular  momentum  in  the  plenum,  which  can  produce  relatively 
large  values  of  angular  velocity  at  the  throat  (e.g.,  in  a  hypersonic 
nozzle  where  r->-0  at  the  throat  on  all  streamlines).  Thus  a 
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fully  justified  approach  to  the  axisymmetric  problem  must  also  include  the 
angular  momentum  equation  in  the  system.  Currently  available  one-dimensional 
analyses2*3  cannot  account  for  either  swirl  or  radial  nonuniformities  of  inlet 
flow  properties,  and  two-dimensional  analyses  appear  to  be  limited  to  supersonic 
reacting  flows4  or  perfect  gases  5'6. 

In  response  to  these  requirements,  General  Applied  Science  Laboratories,  Inc., 
has  developed  Program  R2D2,  a  Fortran  code  for  two-dimensional  chemically  re¬ 
acting,  vibrational ly  excited,  hyperthermal ,  inviscid  flow  in  a  duct,  which  in¬ 
cludes  swirling  axisymmetric  flow.  The  program  name  R2D2  refers  to  its  evolu¬ 
tion  from  two  pre-existing  codes;  one  for  one-dimensional  steady  reacting  flows 
and  the  other  for  two-dimensional,  unsteady  swirling  duct  flows.  Two  versions 
of  the  program  have  been  produced.  The  f i rst, R2D2NE,  is  for  general  nonequili¬ 
brium  flows,  which  may  include  equilibrated  vibration  as  well  as  nonequilibrium 
vibrational  excitation.  The  second,  R2D2EQ,  is  for  fully  equilibrated  flows. 


SECTION  2 


GENERAL  OUTLINE  OF  THE  APPROACH 


1.  TIME-DEPENDENT  METHOD 

The  technique  of  determining  the  steady  state  solution  of  a  mixed  flow  (sub¬ 
sonic,  transonic  and  supersonic  regions)  as  the  asymptote  of  an  unsteady  process 
is  now  well  established.  The  present  approach  can  be  considered  an  extension  of 
the  one-dimensional,  chemically  reacting  flow  analysis  of  Anderson3,  for  example, 
to  two-dimensional  flows.  The  widely  used  explicit,  second-order  accurate,  finite 
difference  algorithm  developed  by  McCormack7,  also  used  by  Anderson3,  has  been 
adopted  in  view  of  its  simplicity  and  reliability.  Although  more  efficient  impli¬ 
cit  algorithms  are  gaining  popularity,  their  application  to  large  systems  of 
coupled  equations,  as  are  encountered  in  multicomponent  reacting  gases,  is  still 
the  subject  of  current  research. 

The  selection  and  application  of  inlet  and  exit  boundary  conditions  for  internal 
flows  warrants  careful  exposition,  particularly  where  the  flow  is  choked  by  a 
throat  in  the  duct.  Consideration  of  a  chemically  reacting  flow  places  an  ad¬ 
ditional  constraint  on  the  proper  boundary  conditions,  the  selected  boundary 
conditions  should,  of  course,  yield  a  unique  as  well  as  steady  solution.  Thus, 
they  must  allow  the  influence  of  intia I  conditions  to  decay  with  the  passage  of 
time,  such  that  the  asymptote  is  indeed  independent  of  the  initial  conditions. 
Finally,  the  sele<  ed  finite  difference  grid  should  conform  to  the  duct  wall  con¬ 
tours,  along  which  the  sole  boundary  condition  of  impermeability  is  to  be  en¬ 
forced.  In  addition,  it  should  afford  a  degree  of  flexibility  in  matching  the 
grid  spacing  to  the  geometric  variations  of  the  duct.  The  present  approach  to 
meeting  these  general  requirements  of  the  time-dependent  method  are  outlined  as 
fol lows . 

2.  BOUNDARY  CONDITIONS 

It  will  be  seen  in  the  following  section  that  casting  the  system  of  equations  in 
characterist ic  form,  at  least  approximately,  provides  useful  guidelines  for  deter¬ 
mining  the  number  and  type  of  boundary  conditions  to  be  imposed  at  the  boundary 
surfaces  of  the  computational  domain.  The  characteristic  form  is  also  often  useful 
in  implementing  a  second-order  accurate  solution  algorithm  at  the  boundary  points. 

The  wall  surface  boundary  condition  is  well  known;  impermeability  requires  a 
vanishing  component  of  velocity  normal  to  the  surface,  the  normal  component 
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of  velocity  must  also  vanish  at  an  axis  of  symmetry.  However,  the  inlet  and 
exit  boundary  conditions  are  not  so  easily  dismissed.  It  ca  be  shown5'8  that  at 
a  subsonic  inflow  boundary  (i.e.,  flow  crossing  into  the  computational  domain) 
all  flow  properties  except  one  must  be  specified  £  priori,  and  at  a  supersonic 
inflow  boundary  all  must  be  specfied.  In  the  former  case,  one  property  must  be 
allowed  to  change  in  time  in  accord  with  the  requirements  of  the  upstream 
travelling  waves  crossing  the  boundary.  In  a  choked  flow,  the  waves  reflected 
off  the  throat  must  be  allowed  to  change  the  inlet  mass  flow  rate  until  a  steady 
state  is  reached.  In  an  unchoked  subsonic  flow  the  waves  reflected  off  the  exit 
boundary  perform  the  same  function.  However  the  thermochemical  state  of  the  in¬ 
let  flow  must  also  be  specified  as  a  boundary  condition.  Therefore,  the  thermo¬ 
dynamic  properties  at  the  inlet,  e.g.,  the  density  or  temperature  cannot  be 
allowed  to  change  in  time.  Consequently,  the  axial  velocity  must  be  selected 
as  the  inlet  flow  property  to  be  determined  by  the  outgoing  waves. 

It  is  also  easily  shown5,8  that  at  a  supersonic  outflow  boundary  (i.e.,  the  flow 
crossing  out  of  the  computational  domain)  no  boundary  condition  can  be  imposed, 
and  at  a  subsonic  outflow  boundary,  only  one  boundary  condition  can  be  imposed. 

In  the  latter  case,  the  static  pressure  is  the  usual  choice. 

In  an  unchoked  subsonic  duct  flow,  this  combination  of  boundary  conditions  does 
not  yield  a  unique  solution,  since  a  variety  of  mass  flow  rates  (i.e.,  axial 
velocities)  can  satisfy  the  stated  conditions.  (This  is  not  an  unusual  inviscid 
flow  situation.)  In  this  case  one  can  expect,  at  best,  the  solution  closest  to 
the  initial  conditions.  However,  in  a  choked  flow,  the  mass  flow  is  fixed  by  the 
throat  size,  and  uniqueness  follows.  In  a  supersonic  flow,  uniqueness  is  again 
guaranteed,  even  if  the  flow  recompresses  to  subsonic  speed  at  the  exit  (e.g., 
the  back  pressure  will  control  the  normal  shock  position). 

3.  INITIAL  CONDITIONS 

The  solution  may  be  initialized  in  one  of  two  ways.  If  no  other  information  is 
available,  a  quas i-one-dimens ional  solution  can  be  generated  by  Program  R2D2 
assuming  uniform  total  pressure,  total  temperature  and  Mach  number  at  each  station 
in  the  duct  and  iterating  the  Mach  number  until  a  specified  initial  mass  flow  rate 
is  matched.  (There  is  no  guarantee  that  it  can  be  matched  to  within  an  arbitrary 
tolerance  in  the  vicinity  of  a  throat.  In  this  case  the  specfied  initial  mass  flow 
rate  must  also  be  iterated  to  match  the  choking  mass  rate  as  closely  as 
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desired.)  The  resulting  velocity  vectors  are  then  rotated  to  satisfy  the  wall 
and/or  axis  boundary  conditions,  and  their  directions  linearly  interpolated  be¬ 
tween  these  boundaries.  The  iteration  procedure  assumes  a  dropping  pressure  in 
the  downstream  direction  to  start  the  iteration  at  each  station.  Thus  it  usually 
converges  on  the  supersonic  branch  downstream  of  a  throat;  however,  the  desired 
branch  cannot  be  guaranteed  since  both  subsonic  and  supersonic  solutions  can 
satisfy  the  mass  flow  criterion.  In  some  cases,  a  small  change  in  ths  specified 
initial  mass  flow  rate  can  change  the  branch  on  which  the  procedure  converges. 

In  addition,  it  must  be  recognized  that  this  procedure  ignores  the  boundary  con¬ 
ditions  to  be  imposed  at  the  inlet  and  exit  stations.  Thus  the  boundary  condi¬ 
tions  are,  in  fact,  imposed  impulsively  on  the  first  time  step.  The  magnitude 
of  the  impulse  should  be  minimized  by  consistent  determination  of  initial  mass 
flow  rate,  total  pressure  and  total  temperature. 

The  initialization  procedure  for  an  axisymmetric  flow  includes  conservation  of 
angular  momentum,  based  on  inlet  boundary  conditions,  along  streamwise  grid 
rows,  which  are  reasonably  good  approximations  to  streamlines.  In  the  non¬ 
equilibrium  version  of  the  program,  the  species  concentrations  and  vibrational 
temperatures  are  initialized  by  conserving  their  inlet  boundary  values  along 
streamwise  grid  rows.  The  validity  of  the  assumption  of  a  frozen  flow  as  the 
initial  state  is  highly  dependent  on  the  particular  duct  geometry  and  inlet 
pressure  and  temperature.  A  rapid  transient  in  the  thermochenical  state  fre¬ 
quently  ensues  due  to  this  aspect  of  the  initialization,  which  could  be  remedied 
(in  the  future)  by  provision  for  an  equilibrium  option  in  the  initialization. 

The  equilibrium  version  of  the  program  initializes  to  a  complete  equilibrium 
state,  based  on  the  initial  temperature  and  pressure  field. 

The  second  procedure  available  for  initializing  the  solution  is  to  use  a 
previously  stored  solution,  which  utilized  the  same  grid  dimensions  and  either 
the  same  or  different  boundary  conditions  and/or  duct  geometry.  This  pro¬ 
cedure  is  intended  primarily  to  permit  restarting  a  solution  which  has  not 
completely  converged  on  a  steady  state.  However,  it  may  be  utilized  to  good 
advantage  in  a  parametric  study,  for  example,  by  initializing  each  solution  to 
a  known  baseline  solution,  or  to  the  preceding  solution  in  a  systematic  varia¬ 
tion  of  parameters. 
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1».  GRID  GENERATION 

The  efficiency  and  accuracy  of  the  finite  difference  solution  Is  clearly  de¬ 
pendent  upon  selection  of  an  appropriate  (and  adequate)  grid  system.  Program 
R2D2  has  adopted  the  basic  grid  generation  scheme  of  the  duct  flow  program 
from  which  It  derives8,  in  large  part  because  it  admits  discontinuities  In 
wall  slope.  It  does,  however,  produce  a  non -orthogonal  grid  which  Is  limited 
to  uniformly  spaced  radial  and  axial  intervals  between  the  boundaries  of  each 
computational  domain(l.e.,  an  upper  wall,  a  lower  wall  or  axis,  an  inlet  sta¬ 
tion  and  an  exit  station).  Since  It  Is  non-orthogonal ,  the  elements  of  the 
Jacobian  of  the  transformation  are  calculated  and  stored  at  each  grid  point. 
Consequently,  the  grid  generation  subroutine  could  easily  be  replaced  by  any 
other  non-orthogonal  (or  orthogonal)  grid  generation  procedure  which  meets  the 
general  requirement  of  equal  spacing  in  a  transformed  coordinate  system  and 
duct  boundaries  mapped  to  the  boundries  of  the  grid  system. 

The  first  step  in  the  present  grid  generation  procedure  Is  to  divide  the  complete 
computational  volume  into  a  sequence  of  contiguous  domains,  within  which  the 
number  of  grid  volumes  and  rows  can  be  independently  specified.  Furthermore, 
the  wall  slopes  need  be  continuous  within  a  domain,  but  may  be  discontinuous  at 
the  interface  between  domains.  The  general  arrangement  of  a  sequence  of  do¬ 
mains  Is  sketched  In  Figure  (1). 


The  number  of  grid  columns  In  each  domain  Is  arbitrary  subject  to  the  limitation 
that  the  sum  of  the  number  of  domains  and  the  total  number  of  grid  columns  in 
all  domains  is  less  than  an  upper  bound: 


I 

•dom  +  T. 


dom 


i  =  l 


"'max:  —  ^bound 


(1) 


where  JboUnd  *s  50  In  the  current  version  of  R2D2. 


The  number  of  grid  rows  In  each  domain  Is  also  arbitrary,  subject  to  the  limit: 

Kmaxj  i.  Kbound  (2) 

where  Kbound  Is  21  In  the  current  ve’rslon  of  R2D2. 


The  physical  Cartesian  or  cylindrical  coordinates,  (x,y)  or  (z,r),  are  trans¬ 
formed  In  an  (£,n)  system  in  each  domain  by  defining 
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FIGURE  1.  SCHEMATIC  OF  A  CONVERGING-DIVERGING  DUCT,  SHOWING  THREE  COMPUTATIONAL  DOMAINS,  WITH  A 
k  x  5  GRID  IN  THE  FIRST,  A  5  x  6  IN  THE  SECOND,  AND  A  6  x  5  IN  THE  THIRD,  AND  A  DIS¬ 
CONTINUITY  IN  WALL  SLOPE  BETWEEN  THE  SECOND  AND  THIRD  DOMAINS. 


xL(v> 

or 

zl  = 

zL(r) 

(5a) 

x*(y) 

or 

ZR  = 

zR(r) 

(5b) 

yH(x) 

or 

rH  = 

rH(z) 

(5c) 

yT(x) 

or 

rT  " 

rT(z) 

(5d) 

are  the  coordinates  of  the  left  and  right  side  (inlet  and  exit)  boundaries 
of  the  i t*1  domain  and  of  the  hub  (axis)  or  lower  wall  and  top  (upper)  wall  of 
the  duct.  The  left  and  right  boundaries  are  assumed  to  be  straight  but  not 
necessarily  vertical.  The  lower  and  upper  walls  of  the  duct  are  specified 
by  a  series  of  arbitrarily  spaced  coordinates,  which  are  fit  by  a  cubic 
spline  curve9  within  each  domain. 


A  rectangular  grid  in  the  transformed  (C,n)  system  is  defined  by: 


AC  j 

'/(Jmaxj 

-1) 

(6a) 

An. 

i 

=  1/(Kmaxi 

-1) 

(6b) 

Obviously,  the  last  grid  line  of  the  ith  domain,  Cj  =  1,  is,  by  definition, 

coincident  with  the  first  line  of  the  next  domain  Cj+i  =  0.  Patching  of  so¬ 

lutions  along  the  interface  between  domains  is  accomplished  by  defining  a  verti 
cal  grid  column  at  C  j  =  -  Af,  .  for  every  domain  (except  the  first),  which  there¬ 
fore  overlaps  into  the  preceding  domain,  as  shown  in  Figure  (2).  Grid  columns 
at  f.  =  0,  Af, . ,  2AC.,  ...  1-AC.  are  considered  "interior  columns",  while  the 

columns  at  C.  =  -  AC .  and  C.  =  1  are  considered  "virtual  columns".  The  solu- 

i  i  i 

tion  at  the  "interior  columns"  is  advanced  by  the  finite  difference  algorithm. 
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<  <J 


p  INTERNAL  GRID  POINT  OP  n™  DOMAIN 
■  EXTERNAL  GRID  POINT  OF  n™  DOMAIN 
o  INTERNAL  GRID  POINT  OF  (n+0™  DOMAIN 
•  EXTERNAL  GRID  POINT  OF  (n+D™  DOMAIN 


FIGURE  2.  GRID  POINT  ARRANGEMENT  AT  INTERFACE  BETWEEN 
n™  AND  (n+D™  DOMAINS 
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while  the  solution  at  the  "virtual  columns"  is  interpolated  from  the  "interior 
columns". 

Note,  first  of  all,  that  the  complete  computational  volume  is  covered  by 
"interior  columns"  of  grid  points;  the  virtual  columns  are  artifices  introduced 
for  convenience  in  patching  the  domains  together.  Second,  a  one-dimensional 
linear  interpolation  is  required  at  points  along  £ .  =  1  if  K^x.  j*  Kfnaxj+1 , 
whereas  a  two-dimensional  linear  interpolation  is  required  at  points  along 

C.  =  -  A£.  if  (VXR)./(J">axi-l)  *  (xL-xR).  +  1/(JmaXj-|)  and  ^max;  ^  •Snaxj  +  i ' 

Third,  the  choices  of 'values  of  Jmaxj,  K^iax;*  *dom  an<*  domain  boundaries 
are  almost  entirely  arbitrary.  This  degree  of  flexibility  is  intended  to  allow 
the  user  to  tailor  the  grid  system  such  that  adequate  resolution  can  be  obtained 
where  it  is  needed  without  burdening  other  areas  with  excessive  numbers  of  grid 
points,  as  well  as  to  accommodate  discontinuities  in  wall  slope.  However,  the 
user  should  be  cautioned  that,  as  a  rule  of  thumb,  large  changes  in  grid  size 
across  a  domain  interface,  and  ratios  of  grid  sizes  in  the  two  directions  sub¬ 
stantially  different  than  unity,  are  to  be  avoided.  Abrupt  changes  in  grid 
size  across  a  domain  interface  tend  to  interfere  with  proper  propagation  of 
waves  across  the  interface  by  filtering  out  and/or  reflecting  the  short  wave 
length  components,  which  cannot  be  resolved  by  the  larger  grid  size.  Large 
differences  in  the  grid  sizes  in  the  two  directions  tends  to  accentuate  the 
numerical  diffusion  in  the  direction  of  the  larger  grid  size,  since  the  smaller 
grid  size  will  control  the  time  step.  In  cases  where  large  ratios  of  grid  size 
are  dictated  by  physical  considerations,  e.g.,  within  a  boundary  layer, 

grid-splitting  techniques  can  be  applied  to  the  basic  McCormack  algorithm  to 
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improve  its  accuracy  and  efficiency  .  However,  grid  splitting  has  not  been 
incorporated  into  the  current  version  of  program  R2D2. 


SECTION  3 

DESCRIPTION  OF  THEORY 


1.  GOVERNING  EQUATIONS 

The  system  of  partial  differential  equations  describing  the  continuum  fluid 
mechanics  of  a  mixture  of  chemically  reacting  and/or  vibrational ly  relaxing 
gases  is,  of  course,  well  established,  cf.  References  (10)  and  (II).  Analytical 
solutions  are  limited  to  highly  idealized  situations,  and  recourse  to  numerical 
solutions  is  generally  required  for  practical  problems.  In  the  present  program, 
the  system  of  equations  is  to  be  solved  by  finite  difference  techniques.  Conse¬ 
quently,  very  little  derivation  is  Involved  in  proceeding  from  the  familiar  form 
in  which  the  system  is  usually  stated  to  the  final  form  in  which  the  differen¬ 
tials  are  replaced  by  differences. 

The  present  program  was  constructed  utilizing  portions  of  a  GASL  program  written 
for  two-dimensional  internal  flow  of  a  perfect  gas8  and  portions  of  a  CALSPAN 
program  written  for  one-dimensional  chemically  reacting-vibrational ly  relaxing 
f low12’ 13’ 14.  Therefore,  it  was  necessary  to  first  adopt  a  consistent  system  of 
nondimensional ization  that  would  be  applicable  to  both  programs.  Thereafter, 
development  of  the  two-dimensional  reacting  duct  flow  program  R2D2,  was  largely 
a  matter  of  logic  and  proper  Interfacing  of  the  components.  Description  of  the 
governing  equations  is  thus  arranged  in  corresponding  fashion. 

a.  EQUATIONS  OF  STATE  AND  NONDIMENSIONAL IZATION 
The  governing  system  of  equations  can  be  completely  nondimensional i zed  by  choo¬ 
sing  two  thermodynamic  properties,  a  length,  L,  a  molecular  weight,  MW*  ,  and  the 
universal  gas  constant,  Rc,  as  the  scale  factors.  The  choice  of  thermodynamic 
properties  is  arbitrary,  but  reduction  of  the  equation  of  state,  which  must  be 
used  repeatedly  in  the  program,  to  a  convenient  form  provides  one  rationale.  We 
have  selected  the  pressure,  P«,  and  temperature,  T„,  as  the  two  basic  thermo¬ 
dynamic  properties. 

A  reference  density  can  then  be  obtained  from  the  equation  of  state; 

o.  -  K.  <» 


11 


and  a  reference  velocity  can  be  defined  by 


a  =  (R  T  /MW  =  (p  /p  >*  (8) 

oo  O  00  00  00  00 

Obviously,  this  reference  velocity  differs  from  a  true  speed  of  sound  by  a 
factor  of  the  square  root  of  the  isentropic  exponent,  y. 

The  resulting  system  of  scale  factors  by  which  variables  are  nondlmensional- 
ized  is  summarized  in  Table  I.  In  particular,  note  that  the  nondimensional 
equation  of  state  is: 

p  =  p  T/MW  =  p  T  £  y.  (9) 

and  the  nondimensional  sound  speed  is: 

a  =  (yp/p)2  (10) 

where  y  is  the  appropr iate I y  defined  isentropic  exponent. 


b.  FLUID  MECHANICS 

The  conservation  laws  for  mass,  momentum  and  energy  for  a  two-dimensional  flow  a 
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TABLE  I 


SCALE  FACTORS  FOR  NOND I  MENS  I ONAL I ZAT I  ON 


VARIABLE 

S 

CALE  FACTORS 

Present 

Reference  (8) 

References  (12)-(14) 

t ,  time 

L/a 

00 

L/a 

oo 

None 

p,  pressure 

p„ 

p»a® 

P. 

p,  density 

II 

p 

00 

Poo 

Poo 

T,  temperature 

T 

00 

a2  MW  /R 

OO  OO  Q 

T 

OO 

v,  velocity 

a 

CO 

a 

oo 

aoo 

H,h,  mixture  enthalpy  and 

E&e,  interna  1  energy* 

a2 

CO 

Q> 

8  ^ 

MW  /R  T 

OO  Q  oo 

h|  ie.,  species  enthalpy 

and  internal  energy* 

R  T 

O  00 

None 

R  T 
o  00 

C  ,C  ,  mixture  specific  heat 
p  v  r 

MW  /R 

00  o 

None 

MW  /R 

00  0 

Cpj,Cv.,  species  specific  heat 

R 

o 

None 

R 

o 

MW,  species  molecular  weight 

MW 

00 

None 

MW 

OO 

w.,  mass  rate  of  production 
of  i  m  species 

L/C.a„ 

None 

L/p  a 

00  oo 

'‘Includes  vibrational  energy 
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3pX2u3  .  3pVVi  .  apx2u2u3 
1 


J  (  3t  +  3x,  +  3x. 


+  P  u2u3 


=  0 


(14) 


3pE  3pu^H  3pu2H  pu2H 

3t  +  3Xj  +  3x2  +  J"  x2 


=  0 


(15) 


where: 


J  =  < 


0  for  planar  flow 


1  for  ax i symmetric  flow 


(x^ ,x2 ) 


(  (x,y) 

<  > 


(z,  r) 


for  j  = 


/  ) 

0 


1 


(16) 


(17) 


<U1 *u2,u3) 


(v  ,v  ,0) 
x'  y 


(vz'vr've)l 


for  j  = 


(IS) 


E  =  H  -  p/p 


H  =  h  +  —  (u2  +  u^  +  u|) 


(19) 

(20) 


This  system,  as  stated,  is  not  complete  without  some  constituitive  relation¬ 
ship  providing  h=h(p,p).  Toward  this  end,  the  thermodynamic  definition  of 
static  enthalpy  is  introduced: 
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h 


(21) 


=  I  a.  h./MW1  =  l  y. 


T 

J  Yj  |y  Cp.  dT  +  Ah.°  + 


eV ! (Tyj) 
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S  '  Cpi 

(T) 

i 

(22) 

Y .  =  Y  • 

1  1 

(p  >T,t ,Yj ) 

(23) 

H 

< 

II 

-1 

< 

(p.T,t,(Yj  )) 

(24) 

The  dependence  of  y.  and  Tv.  on  time  and  on  the  other  chemical  constituents 
relates  to  nonequilibrium  flows,  in  which  solution  of  the  species  conservation 
equation  is  required  to  determine  Yj,  and  solution  of  the  vibrational  energy 
equation  provides  the  vibrational  temperature,  Tv..  The  dependence  of  Tv-  on 
Yj  is  introduced  by  vibration-dissociation  coupling  models12 .  In  an  equi¬ 
librium  flow  y.=y.(p,T)  and  TV.=T,  but  the  functional  relationship  is  usually 
very  complex  and  therefore  expressed  graphically,  rather  than  analytically. 

The  functional  dependence  of  species  specific  heat  on  temperature  is  also 
typically  non-ana  1 yt ica 1 ,  except  for  very  simple  molecules.  Consequently, 
Equations  (11),  (15)  and  (19)-(21)  seldom  yield  an  analytical  solution  for 
temperature.  Therefore,  an  analytical  relationship  of  the  form  h=h(p,p)  is 
generally  unattainable  (in  either  nonequilibrium  or  equilibrium  flows)  except 
for  idealized  systems.  The  usual  procedure  is  to  determine  T=T(h,Y.)  by  itera¬ 
tive  solution  of  Equation  (21).  Since  the  iteration  would  have  to  be  carried 
out  at  every  grid  point  on  each  time  step,  it  could  represent  a  substantial 
computational  burden  in  the  present  program.  Therefore,  it  is  considered 
highly  advantageous  to  rewrite  the  energy  equation  in  terms  of  the  temperature 
rather  than  enthalpy  or  internal  energy,  to  avoid  this  iteration.  Derivation 
of  the  temperature  form  of  the  energy  equation  may  proceed  from  the  second  law 
of  thermodynamics: 


De 

Dt 
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0 


where  e  =  h  -  p/p 

and  D/Dt  =  3/3t  +  u]  3/SXj  +  u2  3/3x2 
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(25) 

(26) 
(27) 


Utilizing  the  eauation  of  state  and  the  definition  of  static  enthalpy,  this 
equation  becomes: 


C 


v 


DT 

Dt 


p2  Dt 


Z 


y  • 


Dev. 

Dt 


Z  e. 


Dy 


i 

Dt 


(23) 


This  may  now  be  combined  with  the  Continuity  Equation  to  yield: 
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e.  =  h .  -  R  T 
i  i  o 


CP  =  E  DT;  CPi 

cv  =  Cp  -  Rq  Z  Y i  =  T  Yj  Cv  j 


and  Y^  =  Cp/Cv 


(29) 


(30) 

(3D 

(32) 

(33) 


This  form  of  the  equation  is  applicable  to  nonequilibrium  flows,  wherein 
the  terms  involving  the  LaGrangian  derivatives  of  the  vibrational  energy  and 
species  concentration  can  be  replaced  by  rate  expressions.  In  an  equilibrium 
flow  the  corresponding  equation  can  be  obtained  by  introducing  an  isentropic 
exponent  defined  by 


a  1  n  p  j 

(3*0 

Ye  =  j 

3  In  p I 

I  / S=constant 

or 

P  = 

k  pY<5  exp(S) 

(35) 

where  the  entropy,  S,  satisfies  Equation  (25)  if: 


DS  =  0 

Dt 


It  then  follows  that 
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(36) 


(37) 


Equations  (29)  and  (37)  are  adopted  as  the  nonequilibrium  and  equilibrium  forms 
of  the  Energy  Equation,  respectively.  It  should  be  recognized,  however,  that 
Equations  (29)  and  (37)  are  restricted  to  shock-free  flows,  whereas  the  origi¬ 
nally  stated  form.  Equation  (15)  is  not.  In  addition,  with  respect  to  the  above 
equation,  it  should  be  pointed  out  that  the  mixture  molecular  weight  only 
changes  significantly  if  the  molecular  composition  of  the  mixture  undergoes  rather 
severe  changes,  e.g.,  from  primarily  diatomic  to  monatomic  molecules.  Thus, 
under  a  wide  range  of  conditions  of  practical  interest,  e.g.,  air  at  temperatures 
up  to  about  7000°K  and  pressures  above  10  ^  atmospheres,  the  last  term  in  Equa¬ 
tion  (37)  is  negligible. 


c.  NONEQUILIBRIUM  CHEMICAL  KINETICS  AND  VIBRATIONAL  RELAXATION 

A  nonequilibrium  flow  may  be  characterized  by  a  time  constant  for  chemical  rate 
processes  and/or  vibrational  relaxation  which  is  of  the  same  order  or  slower 
than  that  for  wave  propagation.  This  characterization  is  subsequently  put  on 
a  more  quantitative  basis  in  the  context  of  stability  considerations  for  the 
finite  difference  algorithm.  If  the  flow  is  out  of  equilibrium  with  respect 
to  chemical  composition  and/or  vibrational  state,  the  preceding  system  of  fluid 
mechanical  equations  must  be  supplemented  with  a  system  of  rate  equations 
describing  the  ponservation  of  chemical  species  and  vibrational  energy.  These 
are  written  as: 


3py . 

It" 


3Pu,T j 
3x, 


SpU2  Yi 
3x„ 


+  j  pu2  Y./x2  = 


MV/. 


=  —  l  Q 
U12  J 


V 


(38) 


i  =  1,2,3... 
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3pev  3pu.  ev .  3pu  ev . 

- -J-  +  - - - J-  +  - - - ^  +  j  pu7ev./x7  =  pev  (39) 

3t  3x1  3x2  Z  vj  z  j 

j  =  f+1 , f+2, f+3, . . .g 

where  the  mixture  of  gases  is  composed  of  I  species,  of  which  those  in  the 
range  f  <  i  <_  g  are  assumed  to  be  undergoing  vibrational  relaxation.  The 
vibrational  states  of  the  species  in  this  range  are  characterized  by  a 
Boltzmann  distribution  ?t  a  temperature  Tvj  /  T.  The  vibrational  energy  is 
thus  expressed  in  terms  of  Tv.  as: 

ev j  =  (0v./TVj  )/(exp(eV|./Tv. )  -  1.0)  (40) 

The  species  in  the  ranges  1  <_  i  <_  f  and  g<  i  _<  I  may  also  possess  vibrational 
degrees  of  freedom,  but  they  are  assumed  to  be  fully  equilibrated,  i.e.,  TVj=T. 

The  mixture  of  I  chemical  species  may  be  composed  of  c  elements  (c>J).  Since 
elements  cannot  be  created  or  destroyed  by  the  reactions,  Equations  (38)  must 
contain  a  subset  for  which: 


'V 

3py  3pu  y  3pu?  y 

- L  +  - ! - L  +  - 1 — L 

3t  3x1  3x^ 


,2,3. . .C 


+  jpu2  ‘/x2  “  0 


where: 


T:  =  E  y 


and  the  6j.  coefficient  matrix  is  composed  of  numbers  of  atoms  of  elements 
in  species  j.  It  therefore  follows  that: 

1  .  c  .  I 

£  5; ;  w  =  0  =  T  5;  .  w ,  +  Z  «i  .  w 


j=1  J  J  j=1  J  J  j=c+1 


i  1 ,2,3. . .c 


1 J  J 


or 


(44) 


-1 


i  =  1 ,2,3. . .c 


where  |  <ScC I  denotes  the  square  matrix  composed  of  the  first  c  columns  of  the 


6jj  matrix,  | 5cj  |  represents  the  remainder  of  the  6jj  matrix  (i.e.,  columns 


c+1  to  I),  and  | w j |  denotes  the  array  w.  from  i=c+1  to  I.  Thus  it  is, in  fact. 


only  necessary  to  evaluate  the  chemical  production  terms  for  species  c+1  through 
I  to  determine  the  values  of  w;  for  all  species. 


Note  that  Equations  (38)  and  (39)  can  also  be  written  as 

Dy.  *i 
ChT  =  pMW . 


(45) 


Dev  j 

’  Dt 


(46) 


or,  for  example, 

Y  j (x j ,x2' t )  =  y.  (x^-Ax^  .x^-Ax^t-At) 


t+At 


t 


Ax^  Ax^ 

where  -  =  — —  =  At 

u,  u0 


(47) 


(48) 


The  latter  form  of  the  equations  is  useful  along  known  stree ,i  surfaces,  such 
as  the  walls  and  axis  of  symmetry  which  automatically  satisfy  the  conditions 
AXj/Uj  =  Ax^/u^. 


The  chemical  production  terms  are  evaluated  from  the  Law  of  Mass  Action  and 
a  set  of  specified  reaction  rates  as  applied  to  a  system  of  chemical  reactions 
in  accord  with  well  known  procedures.  The  reader  is  referred  to  References 
(12)— (14)  for  the  speci f ic  details  of  the  procedure  used  in  the  present  program, 
or  Reference  (15)  and  (16)  for  a  general  overview. 
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d.  EQUILIBRIUM  FLOW 

A  flow  In  which  the  characteristic  time  constants  for  all  relevant  chemical  and 
/or  vibrational  rate  processes  are  much  smaller  than  that  for  wave  propagation 
Is  considered  to  be  in  thermochemfcal  equilibrium.  This  condition  is  in  a  sense 
a  singular  perturbation  of  the  nonequl 1 Ibrtum  formulation  since  the  chemical 
composition,  yj,  and  vibrational  energies,  eV|,  should  not  be  considered  to  be 
functions  of  (x^x^.t),  but  rather  of  pressure  and  temperature  (p,T).  In  an 
equilibrium  flow,  the  fluid  mechanical  system  of  equations  uncouples  from  the 
thermochemical  state  except  for  the  constltuitive  relationship  required  to  de¬ 
fine  ye(p,T)  and  MW  (p,T).  (See  Equations  (9)  and  (37)). 


2-  CHARACTERISTIC  FORM  OF  INLET/OUTLET  BOUNDARY  CONDITIONS 

Assume,  for  the  moment,  that  the  flow  within  a  distance  Ax^of  the  inflow,  or 
outflow,  surface  is  either  in  local  equilibrium,  or  locally  frozen.  The 
energy  equation  then  assumes  the  familiar  form  given  by  Equation  (36): 


S  =  constant 


(49  j) 


on 


(49b) 


which  can  be  combined  with  the  equations  of  continuity  and  axial  momentum 
to  obtain: 


on 


d  I  oq  p  +  y 
dt  -  a 


3 

3x2 


(x^U7) 


U]+a 


dt 


(51'3) 


(50b) 


where  y  and  a  are  the  equilibrium  or  frozen  isentropic  exponent  and  sound 
speed,  respectively.  A  planar  surface  parallel  to  the  (x^.t)  plane  and  trans¬ 
lating  in  the  x^  direction  at  a  velocity  u^  is  referred  to  as  the  reference 
plane  in  which  Equations  (49a)  and  (50a)  can  be  integrated  along  the  lines 
defined  by  Equations  (49b)  and  (50b).  The  reference  plane  and  characteristic 
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linos  are  depicted  in  Figure  (3). 

It  should  be  noted  that  the  angular  momentum  equation,  the  species  and  ele¬ 
mental  conservation  equations,  and  the  vibrational  energy  equation  can  also 
be  written  in  the  form: 


j 


Dt 


0 


(51 ' 


°Yj 

15t 


=  o 


(52) 


□Yj 

dT 


w  .  /p 


(53) 


Dev  j 
Dt 


ev  i 


(54) 


on 


(55) 


Thus,  it  is  clear  that  the  entropy,  angular  momentum,  gas  composition  and 
vibrational  energy  are  all  convective  properties  of  the  flow.  This  statement 
may  be  augmented  by  noting  that  the  curl  of  the  momentum  equation  may  be 
combined  with  the  continuity  equation  to  yield  the  vorticity  transport  equation: 


on 


_D_ 

Dt 


P 


(56a) 


(5>l ) 


Therefore,  the  vorticity  (or  more  precisely  the  ratio  of  vorticity  to  density) 
is  also  a  convective  property  of  the  flow.  Consequently,  a  I  I  these  properties 
of  the  flow  or  equivalent  combinations  of  flow  variables,  must  be  specified  as 
boundary  conditions  at  an  inflow  station,  since  their  characteristic  lines 
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FIGURE  3.  REFERENCE-PLANE  CHARACTERISTIC  SYSTEM 
USED  AT  INLET  AND  DISCHARGE  STATIONS 


-  - - ,„,,,,  (Wf W^Wjj*,UlU«.yjaj,,  ,|,ui  ,  . 

.  - 


(which  are  streampaths)  originate  upstream  of  an  inflow  station  (see  Figure  3) 
Conversely,  none  of  these  variables  may  be  specified  as  boundary  conditions  at 
an  outflow  station,  since  their  characteristic  lines  originate  within  the  com¬ 
putational  domain  upstream  of  an  outflow  station. 


The  radial  momentum  and  energy  equation  can  be  combined  to  yield 


3x 


as  3H  .  r 
3^+J  V 


V 


(57) 


where  £  and  n  are  components  of  the  vorticity  vector  uj.  Therefore,  the  time 
rate  of  change  of  u^  at  the  boundary  points  along  an  outflow  boundary  is  de¬ 
termined  by  the  (known)  flow  properties  along  the  boundary.  Thus,  Equation 
(57)  can  be  integrated  by  the  same  finite  difference  algorithm  used  at  the  in¬ 
terior  points.  If  the  vorticity  is  specified  at  an  inflow  boundary,  then 
Equation  (57)  may  be  solved  along  the  inflow  boundary  points.  However,  if  the 
vorticity  is  not  specified  as  a  boundary  condition,  then  Equation  (57)  must  be 
replaced  by  a  specification  of  u ^  or  (i^/Uj)  a1°n9  the  inflow  boundary. 


Finally,  it  is  clear  from  Equation  (50)  that  if  u^<a  (the  axial  component  of 
velocity  is  subsonic),  then  one  of  the  characteristic  lines  defined  by  Equation 
(50b),  namely  that  which  corresponds  to  the  upstream  traveling  wave  surfaces, 
dXj/dt=Uj-a,  originates  downstream  of  the  boundary,  while  the  other,  namely 
that  corresponding  to  downstream  traveling  waves,  dx^/dt=u^+a,  always  originates 
upstream  of  the  boundary,  as  depicted  in  Figure  (3).  On  the  other  hand,  if 
Uj>a  (the  axial  component  of  velocity  is  supersonic),  then  both  wave  surfaces 
originate  upstream  of  the  boundary  surface.  Since  the  characteristics  which 
originate  outside  the  computational  plane  must  be  replaced  by  corresponding 
boundary  conditions,  closure  of  the  inflow/outflow  boundary  conditon  statement 
depends  on  the  axial  component  of  the  local  Mach  number.  The  required  bound¬ 
ary  conditions  are  summarized  in  Table  2.  The  normally  anticipated  situation 
for  an  arc  tunnel  or  similar  device  consists  of  cases  (la),  a  subsonic  inflow 
boundary  at  the  plenum  entrance,  and  (2b),  a  supersonic  outflow  boundary  at  the 
nozzle  exit.  However,  the  existing  formulation  includes  all  four  (A)  possi¬ 
bilities. 


23 


'At 


TABLE  2.  BOUNDARY  CONDITIONS  AT  THE  INFLOW  AND  OUTFLOW  STATIONS. 


Equation  (50a)  is  integrated  along  the  characteristic  lines  (50b)  which  originate 
within  the  computational  domain,  e.g.,  along  dx^dt^u-a  at  the  inflow  boundary 
and  along  both  dXj/dt=u^  +a  at  the  outflow  station,  to  complete  the  solution  at 
these  boundary  points.  It  is  emphasized  that  this  procedure  avoids  the  need  to 
evaluate  derivatives  at  the  boundary  points  by  non-centered  approximations,  or 
to  devise  heuristic  procedures  such  as  extrapolation  of  interior  point  values 
to  the  boundary  points.  It  does,  however,  assume  constant  entropy  along  stream¬ 
lines  over  a  zone  adjacent  to  the  boundaries  of  thickness  less  than  a  single 
grid  spacing. 

3.  SOLID  WALL  BOUNDARY  CONDITIONS 

The  boundary  condition  of  Impermeability  is  enforced  at  the  walls,  which  re¬ 
quires  the  flow  to  be  tangent  to  these  surfaces.  The  reference-plane  method- 

of-characterl sties  procedure  can  be  applied  along  the  walls  by  reorientation 

(8b) 

of  the  reference  planev  .  However,  for  the  present  problem  the  alternate 
approach  of  using  a  modified  version  of  the 

?.l* 


interior  point  solution  algorithm  has  been  selected?®  Enforcement  of  the  bound¬ 
ary  condition  and  determination  of  the  flow  solution  along  the  walls  is  carried 
out  in  curvilinear  surface-oriented  coordinate  sytem,  in  which  (s,n)  are  co¬ 
ordinates  parallel  and  normal  to  the  surface  ,(vs,vn)  are  the  corresponding  velocity 
components,  and  $  is  the  surface  angle  measured  with  respect  to  the  x^  direction. 
The  transformation  between  coordinate  systems  is  given  by: 


at  the  wal 1 


u,  =  v  cos*  -  v  sin* 
1  s  n  r 


u-  =  v  sin<j>  +  v  cos* 


=  sin*  ±  *  cos*  ± 


The  normal  momentum  equation  is  replaced  by  the  condition  =  0.  The  members 
of  the  system  of  equations  are  solved  by  the  interior  point  algorithm  applied 
in  the  streamline  coordinate  system,  as  explained  in  Section  if. 


4.  AXIS  BOUNDARY  CONDITIONS 

At  the  axis,  the  limit  for  x^-K),  namely: 


is  applied  to  the  interior  point  equations  of  motion  to  derive  the  equations  at 
the  axis  (or  plane)  of  symmetry.  The  axis  equations  are  then  solved  by  the  in¬ 
terior  point  algorithm,  noting  that  at  an  axis  (or  plane)  of  symmetry. 


I  ^  r  j 

—  .  0  ,  ((.; 


and  therefore  as  Xj-^O: 

f(-x2)  -  f(x2),  f  ¥  u2 
u2(-x2)  *  “  u2(x2) 


(64c) 

(65a) 

(65b) 


26 


SECTION  4 


SOLUTION  ALGORITHMS 

1.  INTERIOR  POINTS  AND  STABILITY  CRITERIA 

The  MacCormack  scheme7  has  been  employed  as  the  basic  Interior  point  solution 
method  to  solve  the  equation  system  (II),  (12),  (13),  (14),  (29)  (or  37),  (38) 
and  (39) •  This  method  is  now  well  known  and  widely  used  due  to  Its  combination 
of  second-order  accuracy  and  relative  ease  of  coding,  however,  the  basic  ele¬ 
ments  are  briefly  reviewed  herein. 


The  stated  system  of  governing  equations  are  written  in  the  form: 


3e  3f 
3t  3x 


+  rr2-  +  h=0  =  +A  +B  U-  +C  |S-  +D  I?-  + 


3x 


3t 


3n 


K 


3n 


3C 


(66) 


where  A,B,C  and  0  are  the  components  of  the  Jacobian  of  the  transforma¬ 
tion  from  (x1 ,x2)  to  (£,n),  and  the  arrays  e,f,g  and  h  are: 


e  =  p  (  1,  u, ,  u2,  jx2u3,  T,  y. ,  ey  ] 


(67) 


2 

f=  [p  Uj ,  (pu |  +  p),  pUjU2>  jpu^Uj  x2>pu^T, 


00,7,,  pu,e  ] 
I 


(68) 


g  =  [puju2,(pu2  +  P),jpu2u3x2,  pu2T, 


pVi*  pU2%,  ] 


(69) 


2  2 

h  -  P  (  ju2/x2,  j u j u2/ x2 ,  j(u£  -  u3  )/x2,  ju2u y 


3u.  3u„ 

(T  (V°  (3^  +  +  jWX2 
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+  (EYj^vj  +  ZejiDj/p)  /Cv) 

(jUjY j/^2  "  Wj/pMWj),  Ou2ev/x2  "  *v.)]  (7°) 


(In  the  equilibrium  version  of  the  program,  the  fifth  component  of  the 
h  array  is  redefined  as  3u.  .  3u, 

T<V>  ir  sf)*  h'W*r 

In  addition  v  Is  redefined  as  y. ,  and  the  seventh  component  of  the 
i  1 

above  stated  system  is  dropped.) 


A  "predictor"  value  of  e (x? ,x2>  t+At)  is  obtained  from  a  forward  differ¬ 
ence  in  time  and  non-centered  differences  in  £  and  n-  Using  i,j  sub¬ 
script  to  denote  an  (  C.n  )  grid  point,  and  superscript  n  to  denote  time,t, 
and  e  to  denote  a  "predictor"  value,  the  first  step  of  the  two-step 
algorithm  is: 


n+1 

"V 

e 

>  .j 


i » j 


+  At 
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A(fn  -fn 

M  i,j+1  i, 


,  )/An  + 
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i  +  1  ,J 


.-f° 
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n 

. . . -q. 


n 
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• )  /A£ 


.  )/An  +  d (g? . , 
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The  second  step  consists  of  essentially  averaging  the  time  derivatives 
evaluated  at  t  and  at  t+At  giving  a  "corrector"  value: 


n+1 

f  n 

n+1  . 

m  -m+1 

<bn+1 

e 

i,j 

-*{eu 

+  e  +  At  ( 

i  .j 

■A(fk,j  - 

fi,j-1)/An 

n+1 

<\,n+1 

n+1 

%n+1 

+  B(?.  . 

i  ,J 

-  f.  ,  -)/AC 
i-1, J 

-  g  j  ^  j  _, ) /An 

(72) 

.0+1  n+1  n+1 

jf\j  <\,  f\j 

♦0(9U  Vi.j 

This  algorithm  is  stable  when  the  well  known  CFL  condition  Is  met,  without 
recourse  to  artificial  damping  or  related  filtering  mechanisms.  The  CFL 
condition  may  be  expanded  in  the  context  of  the  present  system  of  equations 
to: 


At  =  (At"1  +  At'*  +  At"!.) 

wave  chem  v i b 


-1 


(73) 


where 


At 


wave 


r  x  »  >  .  •— -i 

^  ^  2  /AxZ  +  Ax^ 


(74a) 


A  t 


chem 


=  £(a>j/p)J 


-1 


max 


(74b) 


A  *v  i  b  “  (  t  .  ] 


mm 


(74c) 


The  Landau-Teller  rate  constant  x.  is  used  for  the  vibrational  timecon- 

t 

stant  which  neglects  the  possible  effects  of  vibrat ion-dissociat ion  coupling. 
These  are  usually  second-order  type  e:fect5,  except  in  extreme  cases,  such 
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as  may  be  encountered  during  an  initial  transient  in  the  solution.  In 
this  case,  vibration-dissociation  coupling  should  be  neglected  until  the 
asymptotic  portion  of  the  unsteady  solution  is  reached. 

Note  that  there  is  no  preferential  order  in  which  the  spatial  differences 
must  be  evaluated  in  the  predictor  and  corrector  steps.  Therefore,  a 
"rotatirg  star"  modification  to  the  basic  algorithm  is  sometimes  employed. ^ 

In  this  modification,  the  order  in  which  points  (i+l,j),  (i-l,j),  ( i ,  j - 1 ) 
and  (l,j+l)  are  used  Is  rotated,  such  that  the  sequence  repeats  every  four  (4) 
time  steps.  While  it  is  true  that  this  procedure  avoids  imparting  any 
preferential  bias  to  the  solution  associated  with  the  order  in  which  deri¬ 
vatives  are  evaluated,  it  is  also  the  experience  of  the  present  investi¬ 
gators  that  it  does  not  measureably  reduce  the  eigenvalues  of  the  stability 
matrix,  and  may,  in  fact,  produce  a  slight  cyclic  behavior  in  the  asymptotic 
solution,  with  a  period  of  four  time  steps.  This  rotating  star  version 
of  the  algorithm  is  incorporated  in  the  program  as  an  option.  Its  utili¬ 
zation  is  left  to  the  judgement  of  the  user. 

2.  INLET  PLANE 

The  inlet  plane  boundary  conditions  have  two  options  fo  the  user  to  choose. 
Either  the  static  pressure  distribution  pin|et(x2)  or  the  total  pressure  dis¬ 
tribution  p  (x  )  may  be  specified.  If  the  static  pressure  is  given,  it  is 

'inlet  2 

accompanied  by  the  static  temperature  distribution  T(x2),  but  note  that  the 
pressure  and  temperature  distributions  should  be  consistent  with  requirements 
of  the  radi,al  momentum  equation.  If  the  total  pressure  P.j.(x2)  's  9'ven  then 
the  total  temperature  field  T-^Xj)  must  be  specified.  A  distribution  of  axial 
velocity,  u^ (k^)  must  also  be  specified,  but  it  is  only  used  if  the  inlet  flow 
is  supersonic  (or  becomes  supersonic).  In  addition,  the  species  concentration 
profile  y, (x.)  and  the  vibrational  temperature  profile  T  (x,)  must  also  be 

12  V  |  2 

given,  as  well  as  the  normal  and  swirl  velocity  components,  u2  and  u^,  or  the 
flow  deflection  angle  u2/u^  ^eu  u2^'  Utilization  of  the  total  pressure 
option  entails  computation  of  the  axial  velocity  Uj(x2,t)  from  the  compati¬ 
bility  equation  by  iterating  on  the  pressure  field  p(x2>t)  at  each  time  step. 
The  density  is  obtained  from: 


30 


(75) 


where 


(PT/RQ  Tt)  MW 


(76) 


Y 


nonequi librium 
equi 1 i bri um 


(77) 


This  option  is  therefore  somewhat  more  time  consuming,  but  presents  an 
advantage  where  wall  curvature  at  the  inlet  station  makes  i t  difficult  to 
estimate  the  static  pressure  distribution  which  satisfies  the  radial  momen¬ 
tum  equation.  On  the  other  hand,  the  resulting  static  pressure  and  temper¬ 
ature  distribution  at  the  inlet  may  not  be  consistent  with  the  specified 
species  concentrations  or  vibrational  temperatures. 

3.  EXIT  PLANE  SOLUTIONS 

If  tbe  outlet  Is  supersonic  (as  In  the  anticipated  applications)  no  boundary 
condition  should  be  Imposed.  The  characteristic  compatibility  equations  give 
all  fluid  mechanical  components  of  the  solution  while  the  thermochemistry 
solution  is  obtained  from  a  La  Grangian  type  method  (tracing  outlet  boundary 
point  streamlines  back  into  the  computational  domain  at  each  time  step). 

For  a  subsonic  exit  condition,  the  exit  pressure  field  must  be  prescribed, 
while  the  rest  of  the  solution  follows  as  above. 


i*.  WALL  POINTS 

As  mentioned  earlier,  the  wall  point  boundary  condition  is  satisfied  by  wri¬ 
ting  the  equations  in  a  streamline  coordinate  system  along  the  wall.  The  nor¬ 
mal  momentum  equation  Is  replaced  by  the  condition  v  “0.  The  normal  derivatives 

n 

In  the  continuity,  streamwlse  momentum  and  energy  equations  involve  only  the 
gradient  of  the  normal  component  of  velocity.  These  terms  are  approximated  by 
a  centered  difference  evaluated  one-half  grid  spacing  off  the  wall.  Second 
order  accuracy  should  be  thereby  retained,  except  in  regions  of  very  small 
radius  of  curvature  (l.e.,  near  points  of  discontinuous  curvature). 
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Since  the  anguiar  momentum,  species  conservation  and  vibrational  energy 
equations  become  exact  integrals  along  a  wall,  i.e.,  see  Equations 
(55). they  are  solved  in  Lagrangian  form: 


rv0  =  constant  «  lrv0](s_v  At)>  (t_At) 


(78) 


y.  =  constant  + 


^  w. /p  dt 


=  ly  +  w.  At/p ] (  At)(r.>r)  (79) 

'  s  » 


At 


e  =  constant  + 
v. 

i 


K 


dt=  [ev.+ev.  At] 


(s-v  At), (t-At) 


(80) 


At 


The  bracketed  quantities  are  evaluated  at  position  v^At  upstream  of  the 
wall  grid  point  in  question  by  linear  interpolation  of  the  known  solution 
at  time  t-At. 


32 


SECTION  5 


CAPABILITIES  AND  LIMITATIONS 


Program  R2D2  is  designed  to  solve  Internal,  inviscid,  chemically  reacting  and 
vibrational ly  relaxing  flows  through  two-dimensional  (planar  or  axisymmetric) 
converging-diverging  ducts.  The  thermochemistry  components  of  the  program  have 
been  derived  from  well  established  theory  and  utilize  subroutines  from  an  exist¬ 
ing  one-dimensional  steady  flow  program12”14.  This,  the  thermochemical  aspects 
of  the  solution  should  possess  a  degree  of  accuracy  consistent  with  the  limita¬ 
tions  of  the  chemical  kinetic  data  and  thermodynamic  properties  and  with  the 
accuracy  of  the  fluid-mechanical  properties  of  the  solution. 

In  regard  to  the  accuracy  of  the  fluid  mechanical  properties  of  the  solution,  it 
should  be  pointed  out  that  few,  if  any,  "benchmark"  type  solutions  are  available 
for  two-dimensional,  subsonic,  swirling,  internal  flows.  Numerical  solutions  are 
seldom  stated  with  sufficient  information  on  boundary  conditions  and  rate  of  con¬ 
vergence  with  respect  to  grid  size  to  permit  comparisons  to  be  made.  Experimental 
data,  on  the  other  hand,  almost  invariably  include  viscous  effects  which  are  not 
incorporated  in  the  analytical  model.  However,  virtually  exact  inviscid  solution 
techniques  are  available17  for  the  interesting,  and  difficult,  case  of  transonic 
flows  near  the  throat  of  nozzles  having  circular  arc  wall  contours,  based  on 
series  expansions  about  sonic  flow.  Comparisons  have  been  carried  out  using  a 
perfect  gas  (i.e.,  frozen  flow)  version  of  the  present  program  for  the  case  of 
an  axisymmetric  duct  having  a  circular  arc  wall  contour  of  radius  R=10r*,  where 
R  is  the  radius  of  curvature  of  the  wall  and  r*  is  the  radius  of  the  throat. 

The  computational  domain  was  defined  to  extend  over  the  axial  range  -0. 19<z/r*<0. 19. 
The  boundary  conditions  imposed  at  the  inlet  station  consisted  of  uniform  total 
Pressure,  uniform  total  temperature,  zero  vorticity  and  zero  swirl.  A  uniform 
back  pressure  corresponding  to  a  Mach  number  of  1.02  was  specified  at  the  exit 
station  to  ensure  that  the  solution  converged  on  the  supersonic  branch  in  the 
diverging  section  (although  in  fact  this  turned  out  to  be  unnecessary  since  the 
solution  never  jumped  to  the  subsonic  branch). 

A  grid  network  having  5  columns  (Az=.0475r*)  and  20  rows  (Ar».05263r*  at  the 
throat)  was  selected  initially.  The  solution  became  asymptotic  after  400  time 


steps  (i.e.,  the  solution  after  500  steps  was  virtually  indistinguishable  from 
that  at  400  steps).  The  Mach  number  distributions  along  the  axis  and  duct  wall 
are  compared  with  the  solution  techniques  from  Reference  (17)  in  Figure  (4). 

Note  that  in  this  case  the  latter  technique  produces  very  nearly  linear  variations 
and  is,  therefore,  highly  accurate(! .e. ,  the  nonlinear  terms  included  in  the  expan¬ 
sion  are  small  but  not  negligible). 

The  differences  between  the  present  solution  and  the  analytical  solution  can  be 
put  in  a  proper  frame  of  reference  by  noting  that  the  maximum  error  in  mass 
flow  rate,  which  occurs  at  the  station  z=.095r*,  is  4/100  of  1%  and  that  the 
average  error  for  the  5  grid  lines  is  1/100  of  1%.  Nevertheless,  the  fact  that 
the  maximum  differences  occurred  at  the  grid  lines  adjacent  to  the  inlet  and  out¬ 
let  stations  motivated  a  second  computation  with  9  grid  columns  and  20  rows,  which 
required  800  time  steps  to  become  asymptotic.  As  can  also  be  seen  in  Figure  (4), 
the  second  solution  was  brought  into  substantially  improved  agreement  with  the 
analytical  solution.  The  maximum  and  average  errors  in  mass  flow  rate  were  both 
halved.  The  radial  distribution  of  Mach  number  is  shown  in  Figure  (5)  at  station 
z=.095r*;  the  agreement  obtained  in  the  second  calculation  is  virtually  perfect. 

Thus  a  high  degree  of  confidence  can  be  placed  in  the  fluid  mechanical  aspects  of 
the  present  technique  if  an  adequate  grid  spacing  is  employed.  However,  the 
judgement  as  to  adequacy  of  the  grid  spacing  must  be  based  on  some  integral 
properties  of  the  solution,  such  as  conservation  of  mass.  It  appears  that  errors 
on  the  order  of  .04%  in  mass  flow  rate  from  station-to-station  may  signal  de¬ 
parture  from  an  "exact"  solution.  The  magnitude  of  the  error  which  is  tolerable 
in  any  particular  case  must  be  determined  in  the  context  of  the  problem. 

The  running  time  required  to  execute  cases  with  a  large  number  of  chemical  species 
may  be  a  significant  deterent  to  the  routine  use  of  Program  R2D2.  The  program  is 
presently  set  up  to  run  a  system  representing  air,  with  three  elements  (oxygen, 
nitrogen  and  argon),  8  chemical  species,  10  chemical  reactions,  and  three  vibra¬ 
tional  ly  excited  molecules.  The  execution  time  on  a  CDC  6600  mainframe  is  on 
the  order  of  0.15  seconds  per  grid  point  per  time  step,  as  compared  to  about 
0.002  seconds  for  a  perfect  gas.  The  amount  of  time  spent  Is  therefore  dis¬ 
proportionately  greater  than  the  increase  in  number  of  equations.  The  calculation 
of  rate  constants,  species  production  terms  and  vibrational  relaxation  terms  is 
believed  to  be  the  primary  contributor  to  the  difference.  Since  these  calculations 
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cannot  be  reduced  substantially  while  retaining  the  desired  generality  in  the 
thermochemical  model,  the  only  savings  can  be  in  reducing  the  number  of  grid 
points  and  number  of  time  steps.  The  former  is  dictated  by  accuracy  requirements 
and  the  latter  by  the  CFL  stability  requirement,  which  couples  it,  in  turn,  to 
the  number  of  grid  points. 


It  is  pointed  out  that  the  time  to  reach  an  asymptotic  steady  state  is  associated 
with  the  time  for  waves  of  both  families  to  travel  between  the  Inlet  and  outlet. 
Thus  a  subsonic  solution  is  generally  the  slowest  to  converge  while  a  supersonic 
solution  (for  which  both  families  of  waves  travel  downstream)  is  faster,  with  the 
required  time  decreasing  inversely  with  Mach  number.  Thus,  the  time  for  a  solu- 
tion  to  a  mixed  flow  problem,  such  as  a  converging-diverging  duct,  to  reach  an 
asymptotic  state  is  controlled  by  the  subsonic  part  of  the  solution.  This  suggests 
that  cases  involving  a  short  plenum  and  a  long,  hypersonic  nozzle,  for  example, 
can  be  most  efficiently  executed  by  treating  only  the  plenum,  throat  and  a  portion 
of  the  nozzle  sufficiently  long  to  ensure  a  supersonic  flow  at  the  exit  plane. 

The  exit  plane  solution  obtained  in  the  asymptotic  limit  for  this  case  can  then 
be  used  as  inlet  station  conditions  for  the  next  segment  of  the  nozzle,  which 
will  then  converge  to  its  asymptote  more  quickly  than  if  it  had  been  treated  as 
part  of  the  first  case.  Alternatively,  a  more  efficient  technique,  such  as  the 
method-of-characteri sties,  can  also  be  used  for  the  hypersonic  expansion. 


The  program  can  accommodate  arbitrary  thermochemical  systems  by  replacing  the 
chemical  kinetic  and  thermodynamic  data  packages,  which  is  discussed  in  greater 
detail  in  the  second  volume  of  this  report.  However,  it  is  specifically  designed 
to  accommodate  a  variety  of  data  packages  representing  frequently  used  thermo¬ 
chemical  systems. 
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SECTION  6 


CONCLUSIONS  AND  RECOMMENDATIONS 


A  program  for  solution  of  two-dimens ional  (planar  or  axisymmetric)  inter¬ 
nal  flows  of  inviscid,  chemically  reacting  and/or  vi brat ional ly  relaxing 
mixtures  of  gases  has  been  developed.  The  program  has  been  specifically 
formulated  for  application  to  arc-driven  wind  tunnels;  it  therefore 
includes  consideration  of  swirl  as  well  as  axial  and  radial  velocity  com¬ 
ponents,  and  of  nonuniform,  subsonic  plenum  conditions.  The  chemical 
system  presently  incorporates  an  8  species  model  of  air,  with  10  chemical 
reactions  and  3  vibrat ional ly  relaxing  molecules.  However,  the  program 
has  been  structured  to  readily  admit  other  chemical  systems  by  redefini¬ 
tion  of  the  thermodynamic  and  chemical  kinetic  data.  In  addition,  arbi¬ 
trary  planar  or axisymmetric  ducts,  including  discontinuous  wall  slopes, 
can  be  considered.  A  simple  nonorthogonal  grid  is  employed  in  the  finite- 
difference  solution  algorithm,  but  other  grid  generation  schemes  are 
readily  admissible,  so  long  as  they  map  the  duct  boundaries  to  boundaries 
of  the  computational  domain  and  define  the  components  of  the  Jacobian 
of  the  transformation  at  every  grid  point. 

The  conservation  laws  for  unsteady  flows  are  solved  subject  to  steady 
boundary  conditions,  with  the  objective  of  determining  the  asymptotic 
steady  state.  Proper  formulation  of  the  inlet  and  exit  station  boundary 
conditions  is  discussed  for  both  subsonic  and  supersonic  flows.  The  pro¬ 
gram  can  treat  arbitrary  combinations  of  conditions.  Uniqueness  can  only 
be  guaranteed  for  a  choked  subsonic  or  a  supersonic  flow,  however  this 
encompasses  the  intended  application  of  the  program.  In  an  unchoked  sub¬ 
sonic  flow  the  program  will  presumably  converge  on  the  solution  closest 
to  the  initial  conditions,  whereas  in  all  other  cases  it  should  be  inde¬ 
pendent  of  the  initial  conditions.  Provisions  to  restart  the  solution 
with  the  same  or  different  boundary  conditions  have  been  made,  with  a  view 
toward  minimizing  the  computer  time  requirement.  An  explicit  finite  differ¬ 
ence  algorithm  has  been  used,  which  therefore  places  a  severe  restriction 
on  the  permissible  time  step.  Consequently,  a  large  number  of  time  steps 
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may  be  required  to  reach  an  asymptotic  solution,  which  can  require  very 

substantial  amounts  of  computer  time.  The  time  requirement  depends,  of 

course,  on  the  grid  size,  the  chemical  system,  the{  boundary  conditions 

and  the  initial  conditions.  However,  execution  times  on  the  order  of 

.15  seconds  per  grid  point  per  time  step  have  been  encountered,  with  the 

considered  air  model  on  a  CDC  6600  mainframe  computer  system.  Since  the 

2 

total  number  of  grid  points  may  be  of  order  10  and  the  total  number  of 
time  steps  of  order  10^,  the  total  running  time  on  this  computer  can  be 
of  the  order  of  four  hours. 

In  view  of  the  substantial  running  time  the  present  version  of  the  program 
requi res, adopt  ion  of  an  implicit  solution  algorithm  is  recommended. 

Removal  of  the  CFL  stability  restriction  associated  with  the  explicit 
method  presently  used  could  probably  permit  an  order  of  magnitude  reduction 
in  the  number  of  time  steps  needed  to  reach  an  asymptotic  state.  On  the 
other  hand,  implicit  techniques  couple  the  system  of  equations,  and  in 
a  chemically  reacting,  vi brat ional ly  relaxing  flow  the  matrix  operations 
may  partially  offset  the  reduction  in  number  of  time  steps. 

The  present  flow  model  is  also  clearly  idealized  with  respect  to  its  in¬ 
tended  application  in  the  sense  that  it  neglects  silch  important  phenomena 
as  turbulence,  liquid  phase  material  from  the  melting  arc,  particulate 
matter  in  the  air  stream,  ionization  nonequi  1  ibriuiii,  etc.,  as  well  as 
asmuthal  nonuniformities  in  the  flow.  These  phenomena  could  also  be  in¬ 
corporated  in  the  program,  although  at  the  expense  of  additional  computer 
resources,  in  varying  degrees.  However,  the  need  to  make  the  present 
version  of  the  program  a  more  practical  tool  for  routine  use  appears  to 
be  the  first  priority.  Nevertheless,  it  should  be  clearly  recognized 
that  this  program  possesses  a  two-dimensional  capability  which  one¬ 
dimensional  analyses  cannot  provide  (especially  with  regard  to  swirl 
effects  and  abrupt  cross-sectional  area  variations),  but  which  inevitably 
requires  more  substantial  computer  resources  than  a  one-d (mens Ional  analysis. 
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